Spatial solitons under competing linear and nonlinear diffractions 
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We introduce a general model which augments the one-dimensional nonlinear Schrodinger (NLS) 
equation by nonlinear-diffraction terms competing with the linear diffraction. The new terms contain 
two irreducible parameters and admit a Hamiltonian representation in a form natural for optical 
media. The equation serves as a model for spatial solitons near the supercollimation point in 
nonlinear photonic crystals. In the framework of this model, a detailed analysis of the fundamental 
solitary waves is reported, including the variational approximation (VA), exact analytical results, 
and systematic numerical computations. The Vakhitov-Kolokolov (VK) criterion is used to precisely 
predict the stability border for the solitons, which is found in an exact analytical form, along with 
the largest total power (norm) that the waves may possess. Past a critical point, collapse effects 
are observed, caused by suitable perturbations. Interactions between two identical parallel solitary 
beams are explored by dint of direct numerical simulations. It is found that in-phase solitons merge 
into robust or collapsing pulsons, depending on the strength of the nonlinear diffraction. 

PACS numbers: 42.65.Tg; 42.70.Nq; 05.45.Yv 



I. INTRODUCTION 



It is well known that bright spatial solitons in optical media are supported by the balance between the self-focusing 
nonlinearity and diffraction. A great deal of attention has also been recently attracted to more special cases, when 
the diffraction is subject to various forms of management^ so that the diffraction coefficient gradually vanishes or 
periodically changes its sign In particular, it was shown in Ref. [2] that related to such a setting with a vanishing 
diffraction coefficient is another recently proposed class of models, in which stable bright solitons are supported by 
the interplay of the ordinary diffraction and self-defocusing nonlinearity, whose local strength grows at r ^ oo at any 
rate faster than r^, where D is the spatial dimension 

A physically important realization of the effectively weak diffraction, which may therefore be sensitive to corrections 
that are usually negligible, is the setting corresponding to the supercollimation in photonic crystals. This setting was 
demonstrated experimentally [4] and studied in detail theoretically [5|-l2|- Similar situations, and the possibility of 
the existence of solitary waves in them, were also analyzed in other photonic media, such as waveguiding arrays [H, 
cavities [9[, and chains with the quadratic nonlinearity [lo|- Furthermore, it was proposed that similar schemes may 
be implemented as well for acoustic waves in sonic crystals [ll], and in Bose-Einstein condensates loaded into optical 
lattices [12]. 

The smallness of the diffraction coefficient near the supercollimation point suggests that the corresponding model 
should be extended by adding nonlinear corrections to the diffraction term in the usual nonlinear Schrodinger (NLS) 
equation. Such a generalization of the NLS equation was proposed in Ref. [7]. In the framework of the generalized 
model, new properties of spatial solitons were predicted — in particular, the existence of a maximum total power, 
above which solitary beams do not exist, as well as inelastic interactions between the beams. The objective of the 
present work is to extend the model proposed in Ref. [7] , and study basic properties of solitary waves in the extended 
system. The main modifications, in comparison with the previously studied model, are that we propose the most 
general nonlinear-diffraction terms which agree with the Hamiltonian structure of the extended NLS equation (the 
term introduced in Ref. cannot be derived from a natural Hamiltonian), and that these general terms depend on 
two irreducible parameters, rather than a single parameter, which was the case in Ref. [7]. 

The extended model is introduced in Section II, where we also present a variational approximation (VA) [13] for 
the study of fundamental solitary waves in this model, and some additional analytical results. Numerical results for 
the shape and stability of the solitons, including the comparison with predictions of the VA and interactions between 
solitons forming out-of-phase and in-phase pairs, are reported in Section HI. The paper is concluded by Section IV. 



2 



II. THE MODEL AND ANALYTICAL FRAMEWORK 



Our aim is to derive a general NLS model including different contributions to the nonlinear diffraction. This is done 
by developing an extended version of the model proposed in Ref. [7] for the description of the supercollimation in 
photonic crystals. We aim to extend the model by using a Lagrangian formulation that produces terms similar to the 
one introduced in Ref. 0, along with additional terms invoked by the underlying Lagrangian/Hamiltonian structure. 

Our starting point is the Lagrangian density for an extended NLS equation including the nonlinear diffraction in 
a general form, £ = (i/2) — qq^) — where q is the local wave amplitude, the asterisk stands for the complex 
conjugation, ^ is the propagation distance, 77 is the transverse coordinate in the planar waveguide, and the Hamiltonian 
density is 
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with real coefficients j3 and 7. Varying the Lagrangian, L ■ 
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Cdrj, gives rise to the generahzed NLS equation, 
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This constitutes an extension of Eq. (3) from Ref. [7]. In particular, the term — (/3/2) q*q^^ which did not appear in 
Ref. (tI, emerges along with its counterpart, — (/3/2) \q\^qnn^ that accounted for the nonlinear diffraction in Ref. 
in the framework of the self-consistent Lagrangian/Hamiltonian formulation. Additionally, the terms proportional 
to 7, which were absent in the equation adopted in Ref. [^, represent contributions to the nonlinear diffraction at 
the same order as the terms ~ /3, hence they should be included too. The remaining terms in Eq. ([2]) correspond 
to the ordinary NLS equation. Dynamical invariants of Eq. (j2j) are the Hamiltonian, H = Hdr], momentum 
M = i q^qdf], and the total power (norm), which we define as 
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It is relevant to mention another equation including the nonlinear diffraction, which was derived as a continuum 
limit of the discrete Salerno model [14] (see also Ref. [15]): 

iq^ + (1 - m\q\^) + 2 (1 - m) \q\^q = 0, (4) 

with parameter m taking values < m < 1. In fact, it is the same equation as the one proposed in Ref. p?], with 
/3 = m/ (1 — m). As shown in Ref. [HI, Eq. dH conserves the norm and Hamiltonian, 
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However, for the purpose of modeling optical media, these norm and Hamiltonian do not seem natural, unlike those 
given by Eqs. (j3j) and ([T]), therefore Eq. ([2j) provides, arguably, a more appropriate model than Eq. (j4]). 

Following the lines of Ref. (7|, we do not include the fourth-order linear diffraction, that would be represented by a 
term ^ qr]r]r]r] in Eq. (|2]), assuming that the equation a ppli es to relatively broad beams, and our objective is to focus 
on effects of nonlinear diffraction. Further, as in Refs. [l^ 0, we focus here on the case of /3 > 0, which implies 
that the nonlinear diffraction competes with its linear counterpart, as in the opposite case a nonlinear enhancement 
of the diffraction does not lead to particularly noteworthy effects. 

A straightforward virial estimate demonstrates that Eq. ([2]) with /3 > gives rise to a supercritical collapse 
(catastrophic self-compression of the wave field) ; the collapse will be arrested if the fourth-order linear diffraction 
is included too. In other words, it is a weak collapse^ so called because only a small part of the total power (norm) of 
the field is involved into the self-compression, the rest being scattered away (the collapse is called strong^ involving an 
essential part of the total power, in the critical case^ which is driven, for instance, by the quintic self-focusing term in 
the one-dimensional NLS equation [HI, UP]). 

Stationary solutions to Eq. (|2]) are sought in the usual form, q{(,^r]) 
satisfying the ordinary differential equation 



= e*^^Q(7^), where Q{r]) is a real function 
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Rescaling the field amplitude and coordinate as 



Q = VkQ, X = Vkr, 



casts Eq. ([7]) into a normalized form, 
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which is more appropriate for the analysis, as it contains a single free parameter, {f3 + 7)fc. Accordingly, the P{k) 
dependence, following from Eqs. (jH]) and (jS)), takes the form of 
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where P = 7r~^/^ IQ{x)] dx is considered as a function of (/3 + 7) /c. Further, it fohows from Eq. (p!Q|) that 
the cutoff value of the propagation constant, k^o^ which is a border of the sohton stabihty region according to the 
VK criterion, dP/dk = 0, is determined by the equation P + 2 (/3 + 7) k^oP' {{P + 7) ^co) = (the prime stands here 
for the derivative). It fohows from here that the cutoff value of the power, Pmax, above which the solitons do not 
exist (see the top right panel in Fig. [2] below), together with the corresponding soliton's peak power and width (that 
represent, respectively, the largest peak power and smallest width that stable solitons may have), depend on the 
nonlinear-diffraction parameters as follows: 
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where numerical found constants are Ck ~ 0.66, Cp ~ 2.35, Ca ^ 1.28, and Cw ^ 2.42. For comparison, it is relevant 
to mention that in the model based on Eq. (jH the cutoff value of the propagation constant. 
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bounds the existence^ rather than stability, of the solitons [14]. Note that Eq. ([T2]) reveals the same scaling of k^o vs. 
P as Eq. (pTj) demonstrates for k^o vs. (/3 + 7). 

Numerical solutions obtained with the help of Eq. (|9]) are produced in the next section. Before turning to that, we 
present the VA for soliton solutions. This approach is well known to produce relevant predictions for many dynamical 
effects exhibited by solitons [13], such the influence of phase modulations on the formation of solitons [16]. To this 
end, we adopt the Gaussian ansatz with amplitude A and width W ^ 



Q(7?) = Aexp [-T]'' / {2W^)] 



(13) 



which is relevant for small values of the propagation constant k [see Fig. ([T]) below]. With the increase of /c, the 
numerically found solution develops a peakon-like shape, hence we expect the VA to fail at larger k. 
Inserting the ansatz (p!3|) into the Lagrangian density corresponding to Eq. ([7j), 
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where may be replaced by the total power ([3]) corresponding to ansatz (p!3|) . which is P = A^W. This substitution 
yields 
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Then, the Euler-Lagrange equations associated with stationary solutions within the framework of the VA are given 

by 

5Leff _ 1 3(/?-7) P 1 . . 

EHminating P here, we arrive at an equation for the width, 

8kW^ + 6[k{p - 7) - 2]W^ - (/3 - 7) = 0. (19) 

Obviously, this equation can be solved analytically. For instance, the physical solution for /3 = 7 is = 3 /{2k). 

In what follows below, the linear stability of the solitary-wave solutions is analyzed by means of numerical methods, 
with the aim to identify parameter regions where the solutions are dynamically stable or unstable. To this end, 
we consider perturbations of the stationary solutions in the form of q{^,r]) = [Q{r]) + e {v{r]) + iw{r]))] e*^^, with an 
infinitesimally small perturbation amplitude e. Upon the substitution of this into Eq. ([2j) and linearization, we obtain 
the following equations: 

- 7 (2QQr^Wr^ - 3Qf^W + Q'^Wrjr] " 2QQr^r^w) - Q'^W, 
-kv -^Vrjrj ^ P {-Q'^Vrjrj " 2Qr}nQv - Q'^V - 2QQrjVn) 

+ 7 {-2QQr^Vr^ - Q%v - Q^Vr^r^ - 2QQ^^v) + 3Q^v, (20) 

which can be written in a symbolic form as = — Lik;, = L2V^ with accordingly defined operators £1,2 • Looking 
for the stability eigenvalues A, so that v{(_^r]) = e^^v{r]), ^{(.^v) = e^^w{r]), we arrive at the linear-stability problem 
in the form of X^v = —L1L2V, and X^w = —L2L1W. The solution of the linear stability problem is reported in the 
following section. 

A well-known necessary stability condition relevant to fundamental solitary waves in models with the self-focusing 
nonlinearity is given by the VK criterion [13- El, dP/dk > 0. It is shown below that this criterion applies to the 
present model, being also sufficient for the stability of the solitons. 



III. NUMERICAL RESULTS 



We now turn to numerical solutions, and their comparison with results of the VA. Fig. [T] presents the numerically 
obtained solutions of Eq. (|9|) as a function of the propagation constant, /c, for /3 = 1 and 7 = 1 in Eqs. ([7j) and 
(O. To obtain these solutions, the Newton's method was used for solving the respective boundary- value problem, 
with vanishing (i.e., homogeneous Dirichlet) boundary conditions and the derivatives calculated by means of the 
centered-finite-difference approximation. Our initial guess for small /c, i.e., small values of the coefficient in front of 
the nonlinear-diffraction terms in Eq. ([9j), was Q{x) = V^sech x, i.e., the stationary solution of the NLS equation, 
and we used the parametric continuation in k (for each subsequent step in /c, the previous solution was used as the 
initial guess). Naturally, it is observed that, as k increases, the nonlinearity plays a more significant role and the 
relevant solution becomes narrower (i.e., with a smaller width) and taller (with a larger amplitude). 

Figure [2] presents basic properties of the numerically obtained solutions, and the comparison with their counterparts 
produced by the VA. Particular characteristics of the soliton solutions that are displayed here are the amplitude, total 
power, and full width at half maximum (FWHM) [the total power is defined according to Eq. (j3j), cf. Eq. (p!Q]) ]. 
It is observed that, as expected, the VA is accurate for small values of /c; however, it does not provide an adequate 
description for sharply-peaked solutions at larger /c, and it does not predict either the existence border for the solitons, 
that corresponds to the the largest possible value of the total power, P = Pmax, nor the stability-cutoff value of the 
propagation constant, k = /Cco, which are evident in the top right panel of Fig. [2l However, as shown above, these 
values can be found in the exact analytical form given by Eq. (pT]) . which is unrelated to the VA. In fact, the wave 
packets suffer collapse at P > Pmax- It is relevant to mention that curves A{k) and P(A:), which are displayed in the 
left and right top panels of Fig. [2l are qualitatively similar to the respective dependences reported in the inset to Fig. 
2(d) and Fig. 2(b) of Ref. [7], for the equation which is tantamount to Eq. (jl]). 

Two squared eigenvalues which, for sufficiently large /c, give rise to instabilities are also shown in Fig. O as a 
function of k. It should be kept in mind that there is always a vanishing eigenvalue in the system due to the phase 
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(gauge) invar iance. At /c > /Cco = 0.33, when the monotonicity of the power of the solution changes, two eigenvalue 
pairs turn positive (i.e., > 0), leading to instability, in accordance with the VK criterion and with analytical results 
(pT]) . Thus, the sharp-peaked solutions found at large k are unstable. 

Predictions of the stability analysis are confirmed by direct numerical simulations of the evolution of perturbed 
solitary waves. In the top left panel of Fig. [3l one observes that, despite adding perturbations, the solution remains 
intact for small k (where the solitary waves are linearly stable). On the contrary, the self- focusing instability, leading 
to the above-mentioned weak collapse^ which involves a small part of the total power at the center of the beam, and 
splits the rest into diverging jets, is evident for a narrower solution with a larger amplitude (i.e., a larger k) in the 
top right panel of Fig. [3l Strictly speaking, the collapse per se is not correctly described by Eq. (j2j), as it does not 
include the fourth-order linear diffraction (which would ultimately prevent it); however, the equation is adequate to 
describe the evolution of the beam towards the collapse. 

As mentioned before there are two unstable eigenvalues when k > 0.33. One of them corresponds to an odd 
eigenmode, and the other one to an even eigenmode with a positive amplitude at the center. Numerical simulations 
show that, for k > 0.33, the stationary soliton perturbed by the odd eigenmode will end up in the weak collapse, 
see the top right panel in Fig. (|3]). On the other hand, with the even perturbation, the dynamics differ for different 
signs of the perturbations. The stationary soliton with the even-eigenmode perturbation added to it will collapse, see 
the middle left panel in Fig. ([3]). On the other hand, under the action of the perturbation with the opposite sign, 
the solitary wave tends to lower its amplitude by emitting small bursts of radiation and ultimately relaxing into a 
breather, as is illustrated by the middle right panel of Fig. ([3]). Finally, the observation of the effects of the different 
parameter variations in the bottom panels of Fig. (|3|) leads to the conclusion that a higher value of 7 induces earlier 
collapse, while the opposite effect is observed with the increase of /3. 

The above simulations started with the input in the form of sech beams. On the other hand, the VA was based on 
the Gaussian ansatz ([13]), therefore it is relevant to test the evolution of solitons generated by the Gaussian input with 
parameters predicted by the VA. A typical example, for parameters belonging to the stability region, is displayed in 
Fig. m One can observe that small (radiation) jets are shed off at the initial stage of the evolution — obviously, with 
the intention to adjust the Gaussian to the exact solitonic shape — followed thereafter by the stable propagation of 
the resulting wave. 

Finally, we consider the interaction of two solitary beams, for the basic cases of in-phase and 7r-out-of-phase soliton 
pairs (i.e., mutually attracting and repelling ones) . The initial in-phase and out-of-phase pairs without a relative 
tilt (i.e., the solitons with zero relative velocity) were built using two well separated copies of a stationary soliton 
produced by the Newton's method for A: = 1, while parameters /3 and 7, which account for the nonlinear-diffraction 
terms in Eq. (j2j), were varied within boundaries of the stability regime. We observe that, in the presence of these 
terms, the repulsive interaction in the out-of-phase configuration does not change dramatically. However, the two 
parameters affect the interaction in different ways. In particular, the increase of /3 leads to a decrease in the strength 
of the repulsion (or attraction) between the solitons, while the increase of 7 produces the opposite effect, see Figs. [5] 
andO 

While the repulsion, leading to the separation of the solitons in Fig. [5l does not reveal any surprising findings, the 
attractive case is far more interesting. The nonlinear-diffraction terms lead to the loss of elasticity of the collisions 
between NLS solitons, and their eventual merger into a single pulsating beam. Furthermore, it is possible to identify 
two different outcomes of the interaction in this case: The merger may lead to a subcritical scenario, i.e., the appearance 
of a non-collapsing robust pulson, that may slowly relax into a stable stationary soliton through weak radiation losses, 
or a supercritical one, which ends up collapsing. These two cases are represented, respectively, by the middle and 
bottom panels in Fig. [6l In the subcritical case, it is observed that the merger into the stable pulson is delayed by 
the terms cx /3, or accelerated by the ones ex 7, but the merger happens, eventually. On the other hand, in the case 
shown in the bottom row of Fig. [6l the supercritical merger, leading to the collapse, occurs at /3 = 0.12,7 = 0, while 
no collapse happens at ^5 = 0,7 = 0.12 In the latter case, more radiation is observed to be emitted in the course of 
the merger, which, in turn, reduces the total power of the emerging pulson, eventually rendering it subcritical. This 
is a remarkable feature, which is uncommon in models that we are aware of. 

IV. CONCLUSION 

The objective of this work was to introduce a generalized model which accounts for phenomena of nonlinear diffrac- 
tion within the one-dimensional NLS equation from the perspective of the self-consistent Lagrangian/Hamiltonian 
formulation. The so obtained model contains two independent parameters, /3 and 7 in Eq. (|2]). The immediate 
motivation for the development of the model is the application to nonlinear photonic crystals in a vicinity of the su- 
percollimation point, as suggested by Ref. In the framework of the proposed general equation, we have developed 
a systematic analysis of the fundamental solitary waves, using the VA (variational approximation), exact analytical 
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results (wherever possible), and numerical methods, which corroborate our analytical observations for individual soli- 
tons and produce new results for their pairwise interactions. In particular, it has been demonstrated that the VK 
stability criterion is valid for the solitons in the present model, and the solitons exist only up to a maximum value 
of the power. The solitary- wave family consists of stable and unstable parts, the border between which, /cco, and the 
corresponding largest value of the total power of the solitons, Pmax, along with the largest peak power and smallest 
width that the stable solitons attain, were found in the exact form, as written in Eq. (fTTj) . Interactions between 
pairs of identical solitons were studied by means of direct simulations. It has been demonstrated that the merger of 
in-phase solitons may lead to the formation of a robust pulson, if the nonlinear diffraction is not too strong; otherwise, 
the emerging beam collapses. On the other hand, the out-of-phase waves always repel each other and separate, with 
the interaction timescales altered by the nonlinear diffraction effects. 

We expect that such a generalized model may hold promise for future explorations, especially if tuned to experi- 
mental data for relevant photonic crystal structures, so that the relative strength of the parameters considered herein 
can be evaluated in realistic settings. On the other hand, it would be particularly relevant to extend the present con- 
siderations to other regimes including the ones pertaining to self-defocusing nonlinearities, or to higher-dimensional 
photonic crystals (possibly with either sign of the nonlinearity ) . Such studies will be deferred to future publications. 
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FIG. 1: (Color online) The top panel shows a sequence of solution profiles (Q) obtained with the increase of k (recall 77 is the 
spatial coordinate). The middle and bottom panels display the profile of Q, and the corresponding eigenvalues of the linear 
stability problem, for k = 0.1 and k = 0.5, respectively (accordingly, the former soliton is stable, while the latter one is not). 
In this figure and below, the results are displayed for /3 = 1,7 = 1, unless otherwise noted. 
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FIG. 2: (Color online) Characteristics of the soliton family with 13=1, 7=1: The left top panel shows the amplitude vs. /c, 
the right top is the power vs. /c, the bottom left is the FWHM vs. /c, and the bottom right panel displays the two squared 
eigenvalues that determine the instability of the solitary waves. The dashed and solid (blue) lines represent, respectively, the 
predictions of the variational approximation and numerical results. The vertical dotted (red) line marks k = 0.33, at which the 
numerically found total power reaches its maximum, the solitons being unstable at k > 0.33, in the precise agreement with the 
VK criterion and with Eq. (pT]) . 
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FIG. 3: (Color online) Contour plots of the evolution of \q\ with an initial perturbation of amplitude 0.001. The top left panel: 
a stable soliton at /c = 0.3, /3 = 1, 7 = 1. The top right panel: an unstable soliton perturbed by the odd eigenmode at 
k = 0.4, /3 = 1, 7 = 1. The middle left panel: the stationary soliton plus the even perturbation for k = 0.4, /3 = 1, 7=1. 
The middle right panel: the unstable soliton minus the even perturbation for k — 0.4, /3 — 1, 7 = 1. The bottom left panel: 
the unstable soliton plus the even perturbation for k = 0.4, /3 = 2, 7 = 0. The bottom right panel: the unstable soliton plus 
the even perturbation for k = 0.4, /3 = 0, 7 = 2. 
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FIG. 4: (Color online) The same as in Fig. O but for the soliton generated by the Gaussian waveform predicted by the VA, at 
/c^O.l, /3 = 7= 1. 
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FIG. 5: (Color online) The repulsive interaction between out-of-phase solitons. The top panel is for /3 = 0,7 = 0, the bottom 
left for /3 = 0.6, 7 = 0, and the bottom right for /3 = 0,7 = 0.6, respectively. 
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FIG. 6: (Color online) The panels display the interaction of two in-phase solitons. The top panel shows periodic elastic collisions 
of the soliton pair in the integrable NLS equation, with /3 = 0, 7 = (in fact, this pattern corresponds to a breather solution of 
the NLS equation [2Q|]). The middle left and right panels are for /3 = 0.08, 7 = and /3 = 0, 7 = 0.08, respectively. The bottom 
left and right panels pertain, respectively, to /3 — 0.12, 7 = 0, and (5 — 0,7 = 0.12. 



